Generalized power expansions in cosmology 



Alejandro S. Jakubi 

Departamento de Fisica, Facultad de Ciencias Exactas y Naturales, Universidad de 
Buenos Aires, Ciudad Universitaria, Pabellon I, 1428 Buenos Aires, Argentina. 



Abstract 

It is given an algorithm to obtain generalized power asymptotic expansions of the 
solutions of the Einstein equations arising for several homogeneous cosmological 
models. This allows to investigate their behavior near the initial singularity or for 
large times. An implementation of this algorithm in the CAS system Maple V 
Release 4 is described and detailed calculations for three equations are shown. 



1 Introduction 



For several homogeneous cosmological models the Einstein equations, com- 
bined with the constitutive equations for the matter source, can be cast into 
nonlinear ordinary differential equations of the form 

N 

D[y{t)] = Y.A iy B 

i=l 





where y is usually a monotonic function of either the scale factor or the Hubble 
rate, t is the universal time, and Bf are real constants, r is the order of the 
ODE and without loosing generality we can assume that Bf > by eventually 
multiplying the equation by suitable powers of y and its derivatives [1] [2] 
[3] [4]. Equation (1) is well defined provided we restrict to y(t) > 0. This 
restriction appears naturally in the cosmological setting because the scale 
factor is positive, or the Hubble rate is positive for an expanding universe. 

When an exact solution of (1) is not available, one is at least interested in 
obtaining some information about it in the form of an expansion for the limits 
t — > + (usually corresponding to the behavior of the solution near the ini- 
tial singularity) and t — ► oo. However solutions to equations of the form (1) 
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frequently do not have power series solutions with integer or even rational ex- 
ponents (Pusieux series) in these limits. Thus we are led to try "generalized" 
power series expansions of the form 



oo 



i/(f)~Ec/* 



(2) 



3=1 



where Cj and rij are real constants and n\ < n 2 < ■ ■ ■ for t — > + and n\ > 
n 2 > • • • for t — > oo. So c\t nx is the leading term, t ni+1 /t ni — > in either limit 
and the set {t Hi } constitutes an asymptotic scale. Inserting this expansion in 
(1) and performing the necessary asymptotic expansions we get 



where the exponents are real and form an ordered set: e\ < e 2 < ■ ■ ■ for 
t — > + and e\ > e 2 > ■ ■ ■ for t — > oo. So, if the equation (1) admits a solution 
with expansion (2), each of the Ct must vanish and this set of equations fix in 
principle the constants Cj,n,j in (2) up to r — 1 free parameters arising from 
the integration constants of the general solution of (1) (for simplicity in the 
notation the arbitrary constant corresponding to time translation freedom is 
fixed to 0). 

Several critical steps in the search for solutions in the form of generalized 
power expansions involve calculations with large number of terms and this 
number grows very fast with the size of the ODE making hand calculation in- 
convenient. However, as many of steps in these calculations have a mechanical 
nature, use of computer algebra systems appear ideally suited. 

The algorithm of used in obtaining these approximated solutions has partial 
similarity in its initial steps to the algorithm used to test whether a given 
ODE satisfies necessary conditions to has Painleve property, namely that its 
solutions are single valued around movable singularities [5]. That algorithm 
checks whether the behavior of the solutions near the singularities is like (2) 
with integer exponents. However in this paper the study of solutions is re- 
stricted to the positive real axis and no consideration will be made about 
their multivaluedness when extended to the complex plane. 
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(3) 



3=1 \ k=l 
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2 The algorithm 



The objective of calculations is to obtain a truncation of (2) to a finite number 
of terms, say M 

M 

l/w(*) = E c i* BJ ( 4 ) 
J'=l 

The method of calculation is iterative, so that constants cm,^m, M > 1, 
when not free are determined by ci, ni, . . . , cm-i,»m-i- For each step of the 
iteration in M, when inserted jjm in (1), two critical points of the calculation 
are: 

(i) Asymptotic expansion (to order M) of the terms with noninteger Bf. 

(ii) Collect all the terms with the same power of t. 
Thus we arrive at an expression of the form 

D[y M (t)] = J2Dit fl (5) 
i=i 

where Di and /; are real constants. In general the sequence of exponents 
/i, f 2 , . . . , /r is not ordered. Thus, for every step the next task is: 

(iii) Sort the {/i}KK fi . 

This sorting operation is perhaps the most involved part of the whole calcu- 
lation. To describe its role and a procedure to do it, only the behavior for 
t — > + will be considered in the following. 

We start with M — 1, that is yi = C\t nx . When it is inserted in the i-th term 
of ( 1 ) we get a term with exponent 

ft = E (m " h ) B " ( 6 ) 

h=0 

and coefficient 

Ei = AiCp^ B ' ; f[(n 1 -h+ B ' = Vic? (7) 

h=i 

If now we take M — 2, crossed terms appear. Let us take any one of them, 
say all factors from the leading term except that of derivative b. Then we get 
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an exponent 



g'i = 9i + (ri2 - ni) B\ > g { 




As ei is the minimum of the {<7i}i<j<jv> an d sa y that this minimum occurs for 
i G /, so Ci = J2i£i Ei is a function of ci and ni, and we find that only the 
leading term of (2) can contribute to the leading term of (3). 

The requirement that the leading term of (3) is not a constant leads to /ij > 
for % G /. If I has only one element the requirement C\ — implies Ci = 0. 
This implies that / must have at least two elements. In this case it is said 
that these terms balance. This constraint usually determines n\. For instance 
if terms i and k balance and /i, ^ we get 



has real roots. On the other hand, if the terms i and k balance and //j = the 
additional constraint X^=o ^ (-^ 

* - b£) = must hold, Ci remains arbitrary 
and ni is determined by the real roots of the equation v { + v k = 0. 

Each pair of real numbers (ci,rii) obtained from this analysis corresponds to 
the leading term of the expansion of a family of solutions. Thus, subsequent 
calculations must be done separately. 

Both the constants A { and B\ in (1) may contain free parameters of the model. 
Thus, considering that some of the B\ are free, we have to ask for the values 
of the parameters that make these ri\ coincide. They are critical values for the 
parameters as they make the behavior of the solutions change. 

Repeating the argument leading to (8) order by order in M it is easy to see 
that Cm and % only appear in Cj for j > M. This also shows that the first 
M terms of (5), once sorted, are equal to the first M terms of expansion (3). 
In particular, the first M — 1 terms are those already found in step M — 1 of 
the iteration. Thus the main tasks at step M are: 

(iv) Find eM, (if possible) and Cm- 

(v) Solve Cm = for either cm or um- 




(9) 



For this n x we get C\ provided 








(10) 
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Generalizing the expressions (6) and (8) we see that the exponents fi are linear 
functions of rii, . . . , n M 

M N r 

/ ; = E«X + EE« (ii) 

j=l i=l h=0 

where the a\ are in turn linear functions of the B\. To order these exponents 
we need to know first when they can be equal. Consider that fi = fk, then we 
have the equation 

M N r 

E («?' - «i) "i + EE (tf* - /^) 5 " = (i 2 ) 

j=l i=l ft=0 

This equation sets a constraint between the {rij} or equivalently defines an 
(M — 1) -dimensional hyperplane in the vector space {ri\, . . . , %) provided the 
vector Aa ifc = (a, 1 - a£, . . . , af 1 - ajf ) ^ 0. 

Only a sector of the M-dimensional vector space is admissible. First of all, 
the inequalities n\ < n<i < ■ ■ ■ < um must be satisfied. Besides, ni, . . . ,um-i 
have been determined in the previous steps of the iteration. As they are usually 
functions of the parameters contained in (1), the variation of these parameters 
within their range further restricts the admissible sector of the M-dimensional 
space. Within it, if the hyperplane exist, tim becomes a linear function of 
ni, ■ ■ ■ ,tim-i, and from all pairs of exponents in {fi}i <l<R the % must be 
found for which both exponents are minimum. In this way we find tu = fi = 
fk, Um and Cm = A + D k . From Cm = we obtain cm- On the other hand, 
if Aai k = 0, cm remains free and Cm = determines um- 

Within the admissible sector, in each subsector delimited by the hyperplanes, 
the fi can be sorted. The minimum, say f k yields tu and provided Dk has a 
common factor of the form c M and the equation Dk = can be solved for a 
real um, this case also yields the pair (cm,um) with cm free. In all cases only 
real solutions are admissible. 



3 Implementation 

The implementation of the algorithm described in the previous section is made 
in Maple V Release 4. At present, this implementation involves use of some 
routines written in the Maple V programming language as well as interactive 
calculations in the worksheet environment. 



5 



For the task of collecting terms with equal power of t, the command collect 
of the current implementation is not suitable as it cannot collect terms with 
nonrational exponents. So, the first task of the implementation was to write 
some procedures to provide this facility. Namely through a command named 
collect2, given as input an expression, a variable (or an expression), one or 
more options and a procedure it yields, depending on the options, 

(a) a form of the expression collected in generic powers of the variable (assumed 
nonnegative throughout), 

(b) a list of the exponents, 

(c) a list of the powers, 

(d) the (collected) coefficients in the form of a table, indexed by the exponents. 

For options (a) and (d) the procedure given as input is used to output the 
coefficients is a useful way (typically using collect ). The package collect2 
has also some useful procedures that can be used independently like expandl 
(like expand but does not expand exponents), and varsubs (like algsubs but 
allows substitution of subexpressions raised to generic powers). 

The rest of the calculations are made for each step of the iteration at the 
worksheet level. It is planned to automatize some of them when some more 
experience is obtained in solving some failures that show the system. One of 
the main difficulties appears in the sorting the list of exponents when a set of 
inequalities is imposed. Maple provides the assume facility by which proper- 
ties about variables, in this case inequalities, are informed to the system so 
that its routines can make use of this information and take decisions when 
handling these variables. The command is is given to check whether a prop- 
erty is true provided that the necessary properties of the involved variables 
hold. However, in some cases, this command fails in deciding the order of two 
expressions, though mathematically the calculation is well defined and a hand 
calculation can demonstrate the inequality in a few lines. So we use proce- 
dures ordl and ord2, that use numerical comparisons, based on values that 
lay inside the interval of interest. These and some other auxiliary procedures 
are currently collected in a package named general. 

In the following sections three examples will be shown to give a feel of how 
this calculations can be done with Maple. 
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4 Example with three terms 



The system of Einstein equations for homogeneous, isotropic cosmological 
models with a variety of matter sources reduce to the ordinary differential 
equation 

y + yy + Py 3 = (13) 



where (3 is a constant. For instance, the problem of a causal viscous fluid with 
the bulk viscosity coefficient ( proportional to p 1 / 2 and y oc H in the truncated 
Extended Irreversible Thermodynamics theory [6]. Also, the behavior near 
the singularity, when the relaxation term is much more important than the 
viscous term in the transport equation, for generic power-law relation ( = 
ap m [8]. For a time decaying cosmological "constant", A ~ —H 3 A [9]. In a 
phenomenological description of the reheating process in terms of an out-of- 
equilibrium mixture of two reacting fluids [11]. 

Equation (13) is also very interesting from the mathematical point of view as 
it appears in the analysis of the Painleve equations [10], and it is the simplest 
case of a class of nonlinear differential equations that possess form invariance 
under nonlocal transformations, so that they can be linearized and its general 
solution can be obtained in parametric form [12] [13]. 

The general solution of equation (13) is 

y(r 1 ) = [Ae x+V + Be^) 1 ' 2 (14) 

At( V ) =[-% (15) 
J y(rj) 



where A± = (1/2) — 1 ± (1 — 8/3) 1 / 2 are the roots of the characteristic poly- 
nomial of the linear equation and A and B are two arbitrary integration con- 
stants. The analysis of the general solution shows that it possesses, for generic 
(3, movable singularities with asymptotic behavior 



y(t) 



a 
At 



£ c n (/3) 7 n Af 



(16) 



n=0 



where r = 4 — a>0is the Kowalevski exponent ([14]), a = a_ for rj — > — oo 
and a = a + for rj — > oo ((3 < 0) and a± = — 2/A±, c = 1 and 7 is an arbitrary 
integration constant. Thus (13) is a second order differential equation of the 
form (1) for which we can show that 2-parameter families of solutions exist 
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with an expansion like (2). Then it is interesting to see how the algorithm 
described above works in this case. 

We begin by writing equation (13) in the Maple worksheet 

> d:=diff (y(t) ,t$2)+y(t)*diff (y(t) ,t)+beta*y (t) "3; 

and then we begin the iteration. 
4-1 Leading term 

For M — 1, inserting the leading term in (13), we get three exponents 

> f :=collect2(subs(y(t)=cl*t~nl , d) ,t , exponents) ; 

/ := [nl-2, 2nl - 1, 3ral] 

We look for the values of nl such that two or more terms balance 

> nle :=balance(f ,nl) ; 

nle := [-1] 

This shows that all three terms balance simultaneously for nl = —1. 

> f :=collect2(subs(y(t)=cl*t~ (-1) ,d) ,t, exponents) ; 

/ := [-3] 

Then e\ = —3, and cl satisfies a quadratic 

> _coeff [-3] ; 

> els :={solve(" ,cl)}minus{0}; 

Id - cl 2 + l3cl 3 



A 1 + V1-8/3 1 1- VI -8/3, 
ClS :={ 2 p '2 p } 

Thus, two families of real solutions with this leading behavior exist provided 
P < 1/8. 
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4-2 Second term 



For M — 2, using C\,ri\ from the previous step, we look for the exponent 
following —3. 

> f :=collect2(subs(y(t)=cl/t+c2*t~n2,d) , t , exponents) ; 

/ := [-3, n2-2, 2n2-l, 3n2] 

These exponents are linear functions of n2, and become equal only for n2 = 
— 1. As we require that n2 > nl — — 1, we can sort them 

> assume (n2> -1) : 

> sort(f , (a,b)-> is(a<b)); 

[-3, n2 - 2, 2n2 - 1, 3n2] 

and find that e 2 = n 2 — 2. Its coefficient 

> _coef f [n2-2] ; 

c2n2 2 - c2n2 + c2 cl n2 - c2 cl +3(3cl 2 c2 

is linear in c2, so that this parameter remains free. Then we find n2 = 3 — 
cl, where the root n2 = —2 is discarded as it arises from time translation 
invariance of one-parameter solutions. On the other hand, as — 1 < n2, this 
implies cl < 4, so that cl_ is admissible for (3 < 1/8, with — 1 < n 2 < 3; while 
cl + is admissible for (3 < 0, with n 2 > 3. 



^.5 TTwrcZ term 

For M = 3, we look for the third exponent 

> f :=collect2(subs(y(t)=cl/t+c2*t~n2+c3*t~n3,d) ,t .exponents) ; 

/ := [-3, n2 - 2, n3 - 2, 2 n2 - 1, 2 n5 - 1, 2 n2 + n5, + 2 n5, 
+ n3 - 1, 3 n^, 3 nS] 

9 



These are linear functions of n3, where we must take into account that n2 < n3 
and — 1 < n 2 (/3) < oo. Let us see first when these exponents become equal. 

> n3e : =balance (f , n3) ; 



3 1 11 

n3e := [-1, n2, -- - - n2, -3-2n2, 1 + 2 n2, - n2 - -, -2 - n2, 

Z Z z z 



l n2 -^2 + 3n2,^n2 



113, 

3'2 + 2 nl! ) 



These in turn become equal when n2 = —1. 
n2e : =balance (n3e , n2) ; 

n2e := [-1] 



Hence we can sort them 

> assume (n2>-l) : 

> sort(n3e, (a,b)->is(a<b) ) ; 



3 1 1 2 11 

[-3 - 2 n2, -2 - n2, n2, -1, - n2 - -, — + - n2, 

1 ' ' 2 2 ' ' 3 3 22 

2 113 

- n2 - -, n2, - + - n2, 2 n2 + 1, 3 n2 + 2] 

O O Z Z 



We look for the third exponent within each interval. 

> 1:=[ n2, l/2+3/2*n2, 2*n2+l , 3*n2+2 , 3*n2+2+l] : 

> ll:=subs(n2=.3,l) : 

for i from 1 to nops(ll)-l do 

u: = (ll[i]+ll[i+l])/2; 

tabla[i] := [ [1 [i] <n3,n3<l [i+1] ] , 

op (3, sort (f, (a,b)->ord2(a,b,n2=.3,n3=u)))] : 

print (tabla[i] ) ; 

od: 

1 3 

[[n2 <n3, n3 <- + - n2], nS - 2] 
1 3 

[[- + - n2 < nS, nS < 2 n2 + 1], nS - 2] 
z z 

[[2 n2 + 1 < n3, n3 <3n2 + 2], 2n2-\\ 
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[[3 n2 + 2 < n3, n3 < 3 n2 + 3], 2n2-\\ 



and we find that there is balance when n3 = 2 n2 + 1 with coefficient 

> expand(subs (n3=2*n2+l ,n2=3-cl , _coef f [n3-2] +_coef f [2*n2-l] ) ) : 

> simplify (" ,{2-cl+beta*cl~2}) ; 

3 pel c2 2 + 36 c3 + 3 c2 2 + (-17 c3 - c2 2 ) cl +2c3cl 2 



Then —1 < n3(/3) < oo and n3 = 2 (n2 + 1) — 1 so that we recover the first 
three terms of (16) with 

c2 2 (3(3cl +3- cl) 
° 3 ~ 36- Wcl +2cl 2 



proportional to c| as expected from (16). Also, as expected, no other solution 
is found. 



5 Equation with a variable exponent 



In order to treat dissipative processes in cosmology which are not close to 
equilibrium a nonlinear phenomenological generalization of the Israel-Stewart 
theory was developed recently [15]. Processes for which this kind of processes 
may have occured are inflation driven by a viscous stress [15,16], and the 
reheating era at the end of inflation [11]. 

In a spacially flat FLRW universe, Einstein's equations together with state and 
transport equations of the fluid give the evolution equation for the Hubble rate 
[15] ' 




H 2 -« (2H + 3 7 # 2 ) - - 1V 2 H 3 = 



(17) 



where 7, a, v, q and k are parameters describing the thermodynamical prop- 
erties of the fluid. Note that q appears in the exponent. 

Let us insert equation (17) in a new worksheet. 

> d: = (l-k~2/v~2-2*k~2/(3*gainma*v~2)*diff (H(t) ,t)/H(t)~2)* 
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(diff (H(t) ,t$2)+3*H(t)*diff (H(t) ,t) + 

( l-2*gamma) /gamma*dif f (H (t ) , t ) ~2/H (t ) +9/4*gamma*H (t ) ~3) + 
(3*gamma*v~2/ (2*alpha)*(l+alpha*k~2/ (gamma*v~2)* 
H(t)"(q-l))*H(t)"(2-q)*(2*diff (H(t) ,t) + 
3*gamma*H(t) ~2) -9/2*gamma*v~2*H(t) "3) : 

Multiplication by a factor H 3 makes all powers of H positive. 

> d2:=map (simplify, expandl (H(t) ~3*d) , power, symbolic) : 

> dview(d2,t,H) ; 



H 3 H" +3H*H' + -2H 2 {H'f + - tf 6 7 - I HJfjL 

7 4 v 2 

9 H*k 2 H' 3 H 2 k 2 (H') 2 2 H 2 k 2 (H') 2 9 H 6 k 2 7 

2 v 2 7 v 2 v 2 4 v 2 

2 H k 2 H' H" 2 k 2 {H'f 4 k 2 (H') 3 ^ H^ 1V 2 H' 

3 7f 2 3 7 2 f 2 3 7f 2 a 

+ -- l^L + 3H 4 k 2 H +-H e 1 k 2 --H e 1 v 2 

2 a 2 2 



5.1 Leading term 



> f :=collect2(subs(H(t)=cl*t~nl ,d2) , t , exponents , 
x->collect (x, [k, cl , v] ) ) ; 

/ := [6 nl , 4 nl - 2, 5 nl - 1, 3 nl - 3, -nJ (-7 + 6 nl - nl q - 1] 



The values of nl that make these terms balance are 
> nle :=balance(f ,nl) ; 

2 13 1. 

nle : = 0, -1, , - , - , — 

-3 + q -2 + q -4 + q q 

These in turn are function of q, so we look whether there is a critical value of 
Q 



> qe:=balance(nle,q) ; 

qe := [1] 

This shows that q — 1 delimits different behaviors of the solutions as it was 
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already known from investigation of equation (17) [17]. To be concise, in the 
following only the case q < 1 will be considered. Further, sorting must be 
done separately for each interval (— oo, 0), (0, 1), (1, 2), (2, 3), (3, 4), (4, oo), so 
we will restrict to the case < q < 1. 

> nlea: =sort (nle , (a,b)->ordl(a,b,q=.5)) ; 

r 1 3 2 1 

nlea := — , -1, , , - , 

g -4 + g -3 + q -2 + g 



Now we look for the leading exponent within each interval. 

l: = [nlea[l]-l,op(nlea) ,nlea[-l] +1] : 

ll:=subs(q=.5,l) : 

for i from 1 to nops(ll)-l do 

u: = (ll[i]+ll[i+l])/2; 

tabla[i] : = [ [1 [i] <nl ,nl<l [i+1] ] , 

op (1, sort (f , (a,b)->ord2(a,b,q=.5,nl=u)))] : 

print (tabla[i] ) ; 

od: 

1 1 

[[-- - 1 < nl , nl < --], -nl (-7 + q)\ 

[[— < nl, nl < -1], -nl (-7 + q)} 
3 

[f-1 < nl, nl < 1, 6 nl - nl q - ll 

LL -4 + g J 1 

3 2 
[[ < nl, nl < 1, 6 nl - nl q - ll 

2 1 

[[ <nl, nl < ], 3 nl - 3] 

-3 + g -2 + g J ' J 

[[ ^ < nl.nl < 0], 3ni -3] 



[[0 < ni, ni < 1], 3ni -3] 



Thus we find that the leading exponent switch at nl = — 1 and nl = 2/ (—3 + 
q), where terms balance. Let us start with nl = — 1. 

> f :=collect2(subs(H(t)=cl*t~ (-1) ,d2) ,t , exponents , 
x->collect (x, [k, cl , v] ) ) : 

> assume (q<l) : 
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> sort(f , (a,b)->is(a<b)) ; 

[-7+q, -6] 

> _coef f [-7+q] ; 



cl 7 + y cl j_ ^ 

a 2 a 



Thus, Hi — —1 and c\ = 2/(37) gives a leading behavior. Let us see the other 
case. 

> f :=collect2(subs(H(t)=cl*t~(2/(-3+q)) ,d2) ,t, exponents, 
x->collect (x, [k, cl , v] ) ) : 

> sort(f , (a,b)->ordl(a,b,q=.5)) ; 

r 9 ~5 + g -7 + g -13 + g 12 
[ -3 + g' -3 + g' -3 + g ' -3 + g J 

> _coeff [-3*(-5+q)/(-3+q)] ; 

16 1 8 1 16 1 3 

6 cl^ 1V 2 | ( Y 7 7^T# + 377^3!^ ~^7 2 (-3 + g) 3)C 
ct (—3 + q) v 2 



Thus ni = 2/(-3 + g) and 



cl = ( r 



9 7 3 w 4 (9-6g + g 2 ) n( 



4 A; 2 a (-7 + 79 - 2) 



-3+9'' 



give the leading behavior of another family of solutions that, for the sake of 
brevity, we will not pursue further in the step M — 2. No solution is found 
when nl is between the balancing values. 



5.2 Second term for n\ — —\ 



We start by expanding the terms with non integer exponents. 

> d3:=subs({H(t)=cl/t+c2*t~n2, 

H(t)~(5-q)=(cl/t)-(5-q)*(l+c2/cl*(5-q)*t~(n2+l)) , 
H(t)"(7-q)=(cl/t)"(7-q)*(l+c2/cl*(7-q)*t"(n2+l))},d2) : 

Then the exponents are 

> f :=collect2(d3, t , exponents ,x->collect (x, [k, cl , c2, v,n2] ) ) ; 
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/ := [-6, -2 + 4 n2, -4 + 2 n2, -7 + g, -3 + 3 n2, 6 n2, 
-5 + g + 2 n2, -1 + 5 nj8, -5 + n2, -6 + g + rajg] 



These are the balancing values of n2 

> n2e :=balance(f ,n2) : 

> assume (q<l) : 

> sort(n2e, (a,b)->is(a<b)) ; 

r 9 _3 g _4 g _5 g _6 g _7 g _ _1 _ g 
I ^ + 9, 2 + 2 , 3 + 3 , 4 + 4 , 5 + 5 , 6 + 6 , i, 2 2 , 



As we require n2 > —1, only two balancing values are admissible n2 = 
and n2 = — | — |. In the case n2 = — g 

> f :=collect2(subs(n2=-q,d3) ,t , exponents , 
x->collect (x, [k, cl , c2,v,q] ) ) : 

> assume (q<l) : 

> sort(f , (a,b)->is(a<b)) ; 

[-7 + g, -6, -5-g, -4-2g, -3-3g, -2-4g, -l-5g, -6g] 



it turns out that e 2 = — 6. From its coefficient 

> _coeff [-6] : 
expand("/cl~3) : 

simplify (subs (cl=2/(3*gamma) , ")) ; 

Av 2 (a + 3 q 7 (9+1) c2 2^ q-6 3 ( ^ 1} 7 (<?+1) c2 2<-«>) 
3 7 2 a 



we find c2 = ^z^y(^) 9 - In the case n2 = -\ - \ 

> f :=collect2(subs(n=l/2+q/2,d3) ,t , exponents , 
x->collect (x, [k,cl,c2,v,q] )) : 

> assume (q<l) : 

> sort(f , (a,b)->is(a<b)) ; 



[ - 7 + ? '-y + 2'- 6 ' "T 



1 c 93 „ o 

-, — 5 — g, g, —4 — 2 g, 

2 2 2 



we find e 2 = — 13/2 + g/2. Its coefficient however only vanishes for c2 = 0. 
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5.3 Solutions 



We collect here the solutions of eq. (17) found thus far, including those for 
q > 1 whose calculation has not been given in detail in this section. We give 
first the solutions with ri\ = — 1. 

Case 2 — q < n 2 < 1 : 

2 1 ^(7-2)^ + 7 

''I = 7= 7, "2 



+ 1 7 (v 7 ^ + l) 



In this case, (1/3) (l/ (l + V^j) < d < 2/3, and (l - y^) / (l + y/2) < n < 
1. We are assuming that c 2 7^ 0. 



Case n 2 = g < 1: 



2 / 2 v ' 



c > = 3? ^ (19) 



Case n 2 = 2 — q < 1 : 
There are three subcases, 
a) 



2 A: 2 _ _8^(M^ /37y 

Cl ~3 7 A; 2 -w 2 ' ° 2 ~ 9 7 ag (2P - 1; 2 ) I2J ( ^ 



with k > v . 
b) 



ci = — ^ (21) 

3 7 1 + y/2v K J 



and c 2 a long expression of the form N/D, with 

iV = 8V2v* (V2v - if (3 7 /2) 9 (v 7 ^ + l) 



9-3 
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D = [2^2 [(q - 1) 7 - 2] v 4 + [7 (2k 2 - l) (q 

+ [-V2 (l + 2k 2 ) (q - 1) + 2^ (4A; 2 
+ [7(1- 2A; 2 ) (g - 1) - 4A; 2 ] u + v 7 ^ (q 



l)+4(l-k 2 )]2v 3 

l)]v 2 
l)k 2 !} 



ci = 



2 1 

37 1 - v 7 ^ 



(22) 



with v < \[2 and C2 a very long expression that is omitted here. 



On the other hand, for rt\ 



2/{q — 3) we get 



97V (9 - Qq + q 2 ) 
4k 2 a (2 + 7 - 75) 



where q < 1. 



6 Anisotropic universe with a scalar field 

Homogeneous anisotropic cosmological models with a self-interacting scalar 
field has been investigated recently to verify the generality of inflationary 
solutions. For this purpose it has been found convenient to cast the metric 
in the semiconformal form [4]. In these coordinates the Bianchi VIo metric 
becomes 



A scalar field with exponential potential V = V^e^ is interesting from the 
physical point of view because it arises in several particle theories, in the 
effective four-dimensional theories induced by Kaluza-Klein theories, including 
various higher-dimensional supergravity [18] and superstring [19] models [20- 
22]. It is also interesting from the mathematical point of view because it allows 
decoupling of the geometric and matter degrees of freedom [23] . Thus the whole 
evolution of the model is obtained from the solution of a third order ODE for 



ds 2 = e m (-dt 2 + dz 2 ) + G(t) (e z dx 2 + e~ z dy 2 ) 



(23) 



G(t) 




(24) 
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where K — — | and m is an arbitrary integration constant. We will investi- 
gate here the singular behavior of this model, corresponding to G(t) — > for 
t — > t . Then we write eq. (24) in a new worksheet. 

> d:=diff (G(t) ,t$2)~2*G(t)-K*diff (G(t) ,t$2)*diff (G(t) ,t)~2 
-diff (G(t) ,t$3)*diff (G(t) ,t)*G(t)+l/2*dif f (G(t) ,t$2)*G(t)~2 
+m~2*diff (G(t) ,t$2) : 

> dview(d,t ,G) ; 

{G"f G-KG" {G'f - G'" G'G+l G" G 2 + m 2 G" 
6.1 Leading term 

Inserting the leading term in (24) we obtain the exponents 

> f :=collect2(subs(G(t)=cl*t~nl ,d) ,t , exponents) ; 

/ := [3ni -4, 3 Til -2, nl -2] 

The values of nl that make these terms balance are 

> nle :=balance(f ,nl) ; 

nle := [0, 1] 

and the leading exponent inside each interval is 

> 11 :=[-inf inity,op(nle) , infinity] : 

for i from 1 to nops(nle) + l do assume (11 [i] <nl ,nl<ll [i+1] ) ; 
tabla[i] := [ [11 [i] <nl ,nl<ll [i+1] ] ,op(l , sort (f , (a,b)->is (a<b) ) )] : 
print (tabla[i] ) ; 
od: 

[[-oo <nl, nl < 0], Snl -4] 
[[0 <nl, nl < 1], Snl -4] 
[[1 < nl, nl < oo], nl - 2] 

Thus the leading behavior changes at n = 1. Let us consider first that nl < 1. 

> _coeff [3*nl-4] ; 
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cl 3 nl 3 - cl 3 nl 2 - K cl 3 nl 4 + K cl 3 nl 3 



We obtain nl = 1/K and cl remains free. This implies 1/K < 1, that is either 
1/2<K<0oyK> 1. On the other hand, the coefficient for nl > 1 

> _coeff [nl-2] ; 

m 2 cl nl 2 — m 2 cl nl 

shows that this case can occur only for m — 0. In such a case nl remains 
free. We reject the balancing value nl = because it does not yield a singular 
behavior. We note however that G = cl is an exact solution. The other bal- 
ancing value nl = 1 is quite special as G = C\t is also an exact solution with 
Ci arbitrary. 

6.2 Second term for nl = 1 
We obtain the exponents 

> f :=collect2(subs(G(t)=cl*t+c2*t~n2,d) ,t, exponents) ; 

> assume (Kn2) : 

> sort(f , (a,b)->is(a<b)) ; 

[n2 - 2, 2 n2 - 3, 3 n2 - 4, n2, 2 n2 - 1, 3 n2 - 2] 

and we find that e 2 = n 2 — 2 for n 2 > 1, with a coefficient 

> _coef f [n2-2] ; 

m 2 c2 n2 2 - m 2 c2 n2 - c2 n2 3 cl 2 - K c2 n2 2 cl 2 + K c2 n2 cl 2 
-2c2n2 cl 2 + 3c2 n2 2 cl 2 

1 2 

Thus we find n2 = ^+2—K , c2 remains free and the constraint ^+2—K > 

cl cl 

1 must be satisfied. 

6.3 Second term for nl = 1/K 

In this case we obtain the exponents 

> f :=collect2(subs(G(t)=cl*t~(l/K)+c2*t~n2,d) ,t , exponents) ; 
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/:= 
3n2 



.2-2K + n2K l — 2K + 2n2K 



K 

2, 3n2 



K 



, n2 -2, 



l-AK + 2n2K 



4, 



2-4K 



K 



n2K -3 + 2K -1 

' K ' 



K 

2a\ 



K 



The values of n2 that make them balance 
> n2e : =balance (f , n2) ; 



n2e :-- 
3 

K 1 ' 



[II _ 
[ 3 a" 
-1 + 2K 

_ K ' 



K-l , 11 + 2A" 12A" + 3 -1 + 2K 1 
K ' ' ' 3 K ' 3 a" ' K ' A 7 ' 



1 1 



2K 



K 



depend on K, and in turn the values of K that make them balance are 

> Ke : =balance (n2e ,K) : 

> sort(Ke, (a,b)->is(a<b)) ; 

-3 -2-1-111212 3 

' 2 ' ' 3 ' 2 ' 3 ' 4' 3' 5' 2' 3' ' 2' ' ' J 



Of these, only the values in [—1/2, 0) or (1, oo) are allowed. Now we could sort 
the balancing values of n2 within each allowed interval of K where it must be 
taken into account that n2 > 1/K. Clearly this analysis is quite branched and 
we will not pursue it here. We just consider the case of the balancing value 
n2 = 1. 

> f :=collect2(subs(G(t)=cl*t~(l/K)+c2*t,d) ,t , exponents) ; 

-2 + K 1_ -3 + 2K -1 + 2K 

J — I a~ ' a" K K J 



In this case the exponents depend on K with balancing values 
> balance (f ,K) ; 

hi, 1] 



Now we sort the exponents within each allowed interval 

> assume (-1/2<K,K<0) : 

> sort(f , (a,b)->is(a<b)) 
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-3 + 2/T -2 + K -1 + 2K 1 



K K ' K ' K s 



> assume (K>1) : 

> sort(f , (a,b)->is(a<b)) ; 



-1 + 2K -3 + 2K -2 + K 1 
^ K ' K ' K ' ~K 



For K > 1 the coefficient is 

> _coeff [-(-1+2*K)/K] ; 

c2 2 cl m 2 cl m 2 cl cl c2 2 cl c2 2 9 

+ 77- + 3 — - 3 — — + c2 2 cl 



K 3 K 2 K K 2 K 



and we get a solution: n2 = 1 and c2 = ± ^_™ . 



7 Conclusions 



We have shown that approximations of solutions of nonlinear ordinary differ- 
ential equations relevant to cosmology can be obtained with the help of CAS 
in the form of generalized power asymptotic expansions, without much effort. 
For this purpose we have sketched an algorithm to make these calculations by 
an iterative process. 

The implementation of these calculations were made in Maple V Release 4. 
This work environment is quite productive as it allows interactive exploration 
as well as a powerful programming language. The set already written of proce- 
dures to collect terms, obtain exponents and coefficients work in satisfactory 
way. The next step that remains to be coded is sorting of exponents. This 
task is more complex as it involves branching of options and has the further 
difficulty that arises in the weakness of the current version of the assume 
facility. 

From a more theoretical point of view, it deserves to be investigated how the 
complexity of the algorithm increases with the order of iteration, and whether 
this growth puts an effective limit to practical calculations. In such a case it 
would be interesting to know whether more efficient algorithms can be devised. 

Also it would be interesting to know whether expansions of solutions as shown 
in this paper, can give information about the integrability of the equation. 
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